Non-linear oscillatory rheological properties of a generic continuum foam model: 
comparison with experiments and shear-banding predictions 
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The occurence of shear bands in a complex fluid is generally understood as resulting from a 
structural evolution of the material under shear, which leads (from a theoretical perspective) to a 
non-monotonic stationnary flow curve related to the coexistence of different states of the material 
under shear. In this paper we present a scenario for shear-banding in a particular class of complex 
fluids, namely foams and concentrated emulsions, which differs from other scenarii in two impor- 
tant ways. First, the appearance of shear bands is shown to be possible both without any intrinsic 
physical evolution of the material (e.g. via a parameter coupled to the flow such as concentration 
or entanglements) and without any finite critical shear rate below which the flow does not remain 
stationary and homogeneous. Secondly, the appearance of shear bands depends on the initial condi- 
tions, i.e., the preparation of the material. In other words, it is history dependent. This behaviour 
relies on the tensorial character of the underlying model (2D or 3D) and is triggered by an initially 
inhomogeneous strain distribution in the material. The shear rate displays a discontinuity at the 
band boundary, whose amplitude is history dependent and thus depends on the sample preparation. 



PACS numbers: 47. 57. Be Foams and emulsions 83.10.Gr Constitutive relations 
foams in Rhcology - 83. 50. Ax Steady shear flows, viscometric flow 
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I. SHEAR BANDS AND FOAM RHEOLOGY 
A. Shear bands in complex fluids 

It may seem paradoxical that a single material, when 
submitted to a uniform shear stress o~ xy , between two 
parallel plates or two coaxial cylinders, may be observed 
simultaneously in two distinct states in different regions 
of the flow. This 'shear bands' observation has neverthe- 
less become common since the early 1990s in a variety 
of complex fluids: they appear and are stable or 
sometimes fluctuate These bands are most of the 

time parallel with the shearing plates [l[ , with a different 
shear rate in each band. 

The current understanding of these observations relies 
in general on two essential ingredients : (i) a structural 
evolution of the material under shear, and (ii) a stress 
response that decreases as a function of the shear rate 
(within a particular range). This decrease is the mechan- 
ical signature of the structural evolution of the fluid and 



is the source of the mechanical instability that triggers 
the appearance of bands [j| . 

In polymer melts or entangled polymer solutions (lOj 
and in entangled giant micelle solutions @ , the flow elon- 
gates the objects, which alters the apparent viscosity of 
the material (which must be evaluated after subtracting 
the effect of wall slip [H|). The fact that this viscosity 
goes down is principally due to the average orientation of 
the objects in the shear flow. 

In lyotropic lamellar phases, the transition can be asso- 
ciated with the reorganisation of the films into onion-like 
multilamellar vesicle systems [H, 0, [lU, [l~3| , also exhibit- 
ing wall slip behavior [f|. In micellar cubic crystals the 
transition consists in an ordering of the initial polycrystal 
into a single crystal with specific planes becoming aligned 
with the plates [3, [Hi- In the last two cases, no micro- 
scopic interpretation of the decrease in effective viscosity 
occuring during the transition is currently available. 

In granular materials, surface flow is a particular case 
of shear bands: the lower band is in this case blocked 
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(zero shear) while the flowing region is sheared. Again, 
no complete structural description is available. Neverthe- 
less, it is admitted that through dilatancy, which reflects 
the necessity for the grains to move a little bit apart in 
order to move past each other [l6[ , the shear generates a 
difference in volume fraction between the flowing region 
and the blocked one. This lower volume fraction tends 
to facilitate the flow in the flowing region even more as 
compared to the blocked region, thus stabilizing shear- 
banding. When it is present, gravity is of course essential: 
it favours this phenomenon by helping the system segre- 
gate into a dense, blocked region (located at the bottom 
if the particles are denser than the fluid) and a less dense, 
flowing region. Thus, it allows to determine the concen- 
tration profile (l7j . 

In foams and emulsions, the situation is less clear. 
Shear bands were observed in 2D [TH, In some cases, 
the observed shear-banding could result trivially from 
shear stress inhomogeneity, due to cylindical Couette ge- 
ometry (a(r) oc 1/r 2 ) or enhanced (for 2D foams) by the 
presence of at least one solid boundary |2CH23[ whose 
friction on the foam implies that r 2 a(r) is not uniform. 
In some more interesting cases, the shear rate is spatially 
discontinuous at the boundary between the blocked and 
the sheared regions 0, [24| . This discontinuity may arise 
from an apparently intrinsic impossibility for the foam 
to be deformed homogeneously at low shear rates [25| . 
which then implies the presence of shear bands at low 
shear rates. But at least in some cases, the shear rate at 
the boundary is not unique for a given system (2~H and 
is thus not intrinsic. Recently, the very existence of such 
a finite shear rate at the boundary has been seriously 
questioned (2(|. As we shall see, the present work high- 
lights yet another (history dependent) possible origin of 
the shear rate spatial discontinuity, arising from the ten- 
sorial character of the material response. For these ma- 
terials also, no complete structural description accounts 
for flow localization in a satisfactory manner. Dilatancy, 
which corresponds to a local change in water concentra- 
tion (j>, certainly plays an important role by easing the 
relative motion of bubbles or drops, although it behaves 
somewhat differently from granular materials depending 
on the liquid fraction [27H29( | . The structural disorder is 
also invoked as a parameter coupled to the flow [3(J. In 
both cases, the local fluidity (ratio of the shear rate and 
the shear stress) is enhanced. 



B. Foam rheology 



been devoted to address this challenge and describe the 
richer mechanical behaviour of foams. Several rheologi- 
cal models have thus emerged (3ll - [35| . They all assume 
that the foam is essentially incompressible. Within this 
perimeter, some models are purely visco-elastic yet with 
a non-linear elasticity |31| . As for the models incorporat- 



ing plasticity, they can be dispatched into two categories: 
the creep term either depends on the stress and deforma- 
tion rate [H, HH, HU or on the stress only [32|,[34[. Finally, 
these models also differ in the tensorial form of elasticity 
and creep, a feature which is relevant for non strictly 2D 
systems (or for compressible materials). 

Despite this variety of models, most current experi- 
ments are not sufficiently stringent to fully test these 
models and decide which ingredients are indeed relevant 
to describe the mechanical response of foams. For in- 
stance, classical linear rheology experiments, particularly 
oscillatory measurements, are clearly unable to provide 
much information about the behaviour under large stress. 

Because the constitutive objects of foams are macro- 
scopic and can be observed directly, statistical tools have 
been elaborated to measure the local deformation and 
deformation rate. Using these tools, more complex ge- 
ometries such as flows around obstacles are also used in 
order to subject the foam to a tensorially broader variety 
of sollicitations [37|. Yet because these experiments are 



conducted in a confined geometry, the unknown friction 
with the walls and the fact that no direct measurement 
of the total stress is conducted, make it difficult to test 
stress predictions beyond low velocities. 

To complement this, in order to test the full time re- 
sponse of the models even within classical geometries 
such as those available in a rheomctcr, a broad range of 
experiments could be elaborated by choosing many differ- 
ent forms for the time dependence of the applied deforma- 
tion or stress. As a first step towards this goal, consider- 
ing the full time-dependent response of a foam subjected 
to large amplitude oscillatory shear (not only the usually 
extracted storage and loss moduli), is susceptible to pro- 
vide more stringent tests for models. Such experiments 
have been conducted recently: because the flow was ob- 
served to remain homogeneous, the measured behaviour 
can be robustly attributed to the local, yet macroscopic, 
mechanical response of a 3D foam [38| . The strain-stress 
(Lissajous) curves display various shapes and show that 
such data can become available and provide non triv- 
ial results. These results will be discussed later in the 
present work. 



The specificity of foams as compared to other materi- 
als is the following: not only do they flow substantially 
only above some threshold stress, but they also undergo 
large elastic deformations at lower stress. Hence, clas- 
sical models such as visco-elastic fluids (well suited for 
polymeric fluids) and elasto-plastic solids (well suited for 
metals) are unsufficient to capture the behaviour of foams 
and emulsions. In the past few years, much effort has 



C. Scope of the present work 

In materials whose plastic threshold corresponds to a 
small deformation, shear banding requires a fluidizing 
mechanism such as those mentioned in paragraphe II Al 
above. But for materials whose deformation at plastic- 
ity onset is large, like foams, the tensorial nature of the 
material state, due to stored deformation, is sufficient to 



obtain shear bands: no extra mechanism is required. 

In particular, the model that we suggest does not in- 
corporate such an ingredient as dilatancy. We know that 
a local plastic event results in an elastic redistribution 
of stress in the neighbourhood [IH, |4£j . In simple shear 
geometry, it thus favours flow localization [4lM43| . In 
a statistical manner, it then raises locally the material 
fluidity [44j and generates a non-local material rheology. 
This non-local character had been observed in concen- 
trated emulsions flowing in microfluidic channels [45| . 
Let us emphasize that these non-local effects are intrinsi- 
cally present in our modelling since the underlying clastic 
propagators [Zlj - El ] result directly from the elastic con- 
tinuum medium equations that we use. 

The paper is organized as follows. We first discuss sev- 
eral types of continuum models for foams (Section [TTJ) . 
We then recall (Section ITO]) the construction of a rather 
generic continuum model for foam or emulsion rheol- 
ogy 241 • ^ is generic as for the elasticity, the plastic 



flow rate and the specifically three-dimensional form of 
the response. We then simulate large amplitude oscil- 
latory shear (Section IIV[) and conduct a first round of 
comparison with published data on such experiments at 
a fixed frequency for various amplitudes [38j . We are not 
able to reproduce the experimental data in a reasonable 
way with a single set of parameters, but in the future, fit- 
ting similar data over a full range of both frequency and 
amplitude will be a very effective and stringent method 
for testing more general models than the present one. Fi- 
nally, we thoroughly discuss shear banding in our model 
(Section|V|. Note that in this work, we restrain ourselves 
to a strictly mechanical and thermodynamical formula- 
tion. The important problems related to the coupling 
between the rheological behaviour and the structure of 
the material are not discussed. This coupling is experi- 
mentally well documented in various complex fluids sys- 
tems in which shear bands are de facto associated with 
structural transitions . 

Here, we ask a more restricted question: could station- 
ary shear bands in foams and emulsions be accounted for, 
starting from an inhomogeneous initial stress distribu- 
tion in the material with otherwise strictly homogeneous 
mechanical properties? Our main result: shear bands 
can emerge in a structurally homogeneous material un- 
der shear, only due to an inhomogeneous distribution of 
the initial internal stress in the material. We demonstrate 
this for a physically very natural form of the elastic and 
plastic laws. 



II. CHOOSING THE TYPE OF CONTINUUM 
MODEL 

As mentioned above, our main point is to explore mod- 
els without any extra dynamic variables apart from the 
stored local deformation. In the present Section, we will 
discuss these models in very general terms, omitting any 
explicit tensorial features. 
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FIG. 1: Schematic, scalar view of rheological models (defor- 
mation e and stress a) with only one internal degree of free- 
dom, namely the deformation e of the spring, (a) General 
model: it has one (non-linear) spring and three creeping ele- 
ments, each of which may have a flow threshold, (b) Foams 
and emulsions display some (large) deformation before creep 
is triggered. Hence, creep elements 1 and 2 cannot have a fi- 
nite threshold. By contrast, element 3 does have a threshold. 
Here, for simplicity, we assume that the spring and both vis- 
cous elements are linear, and that element 3 cannot withstand 
any stress beyond 03. 



e 






T 


2T\ 


3T 


St 



m/Gi 



Gi 




S'G^ e <^-.. 1 ra 1 




FIG. 2: Response of the model represented in Fig. [Tb to a tri- 
angle wave deformation, (a) Deformation e as a function of 
time, (b) In the limit 772 = 0, stress a (solid line) and elastic 
stress Gi e transmitted by the spring (dashed line) as a func- 
tion of time. When the threshold 03 is reached, the spring 
relaxes with the time scale ni/Gi of the Voigt element, (c) 
Corresponding Lissajous representation of a and G\ e. The 
periodic jump in stress (green segment) has the same ampli- 
tude as the initial jump (red segment) . (d) In the limit 771 = 0, 
stress a (solid line) and elastic stress Gi e transmitted by the 
spring (dashed line) as a function of time, (e) Corresponding 
Lissajous representation of a and Gi e. The amplitude of the 
periodic jump in stress (green segment) is twice that of the 
initial jump (red segment). 
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A. Models with only one internal variable 

The most general arrangement of rheological elements 
having only one internal deformation variable is repre- 
sented in Fig. Q]}. The deformation of the (not necessar- 
ily linear) spring represents the deformation of the local 
structure. Three creep elements (viscous and/or plastic) 
can be included, as shown. 

At this point, let us recall that foams are viscoelas- 
tic under weak stress conditions. Hence, both the spring 
itself and the combination of spring and creep elements 
must be able to deform under weak applied stress. As 
a result, creep elements 1 and 2 must flow under weak 
stress: we cannot choose them with a stress threshold be- 
low which they would not flow at all. In other words, they 
are purely viscous (although not necessarily linear). By 
contrast, creep element 3 must have a stress threshold 
so that the entire system also displays a stress thresh- 
old. These considerations are summarized schematically 
in Fig. [T|d. 

Although elements 1 and 2 are viscous, they play differ- 
ent roles when some creep motion of element 3 is involved. 
Let us first illustrate this point by considering an exper- 
iment in which we impose a constant deformation rate 
from t = and reverse the deformation rate as of t = T. 
For simplicity, we consider a linear spring with modulus 
G\, Newtonian viscosities rji and 772. Furthermore, we 
consider that element 3 is simply a solid friction element 
with threshold 03 with no dependence on velocity. Fig. [3] 
shows the contributions of viscous elements 1 and 2 sepa- 
rately in the response of such a system to a triangle wave 
deformation. In both cases, the stress jumps up immedi- 
ately to a finite value at t = due to the viscous elements 
rji and 772. The stress then increases at a constant rate 
as the spring elongates. When the threshold of element 3 
is reached, the stress saturates and remains constant. At 
t = T, when the deformation rate is reversed, the stress 
jumps down by a finite amount. It then decreases at a 
constant rate as the spring is relaxed and later stretched 
in the reverse direction. 

There are two differences between the effects of viscous 
elements 1 and 2. The first difference is that the observed 
threshold depends on the deformation rate in the case 
of viscous element 2. But that feature is not essential: 
one can always decide that the deformation rate of creep 
element 3 affects its stress (in other words, by considering 
that it contains not only a solid friction element, but also 
an extra viscous element in parallel with it). 

The second difference between both situations of Fig. [2] 
is more essential. The jumps in stress at t = and at t = 
T have equal magnitudes in the case of viscous element 1. 
By contrast, in the case of element 2, the magnitude of 
the second jump is twice as large as that of the initial 
jump (except if T is too short for the system to be able 
to relax the Voigt element, with timescale 771 /Gi, after 
it has reached the threshold stress). More generally, in 
such an experiment, one can express each viscosity (at 
the applied deformation rate) in terms of the magnitudes 
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FIG. 3: (a) Burger model: schematic diagramme (left) and 
deformation response e(t) to a step imposed stress cr(t). 
Spring G'i elongates immediately (red segment) while spring 
Gi responds with some delay (green curve) due to viscous ele- 
ment 771 . Viscous element 773 gives rise to a constant additional 
deformation rate (blue segment), (b) Combined Bingham- 
Burger model: solid friction element 0-3 (complemented by 
viscous element 773) is now included so as to provide addi- 
tional creep above the stress threshold. We believe that some 
tensorial version of this model could mimic the rheological be- 
haviour of a foam quite adequately, (c) Model studied in the 
present work. Apart from the additional viscosity 772 intro- 
duced in Fig. [T] it represents the combined Bingham-Burger 
model when the (blue) viscous element has not been deformed 
yet, either at intermediate time scales when the (green) Voigt 
element has relaxed (the green and red springs then respond 
in series) or at short time scales when the Voigt element is 
still blocked (in which case only the red spring responds). 



of the stress jumps at t = 0, when the deformation rate 
changes from zero to +7, and at t = T, when it is reversed 
from +7 to —7: 



77i(7) 
^2(7) 



2|Aq(0)|-|Aq(r)| 

7 

|Acx(r)|-|A<r(0)| 
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(1) 

(2) 



When the value of the elastic deformation e can be mea- 
sured independently, for instance through optical mea- 
surements in 2D foams and relevant statistical tools 37[ , 
the respective contributions from elements 1 and 2 can 
be obtained by comparing a and e (see full lines versus 
dashed lines in Fig. [2]) . 



B. Weak stress: role of one extra internal variable 

Below the flow threshold, the model outlined above 
behaves like a Voigt element (a spring in parallel with a 
viscous element). The behaviour of a foam under a weak 
stress is in fact a little more complex: a linear Burger 
model (see Fig. [3^) is known to correctly reproduce step- 
wise creep experiments on liquid foams [461 ] . Because 
the Burger model contains two springs, it corresponds 
to a system with two internal variables (spring elonga- 
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tions e and e') rather than one. Fig. [3^ shows the re- 
sponse of such a model to a stepwise creep experiment. 
The stress jump generates an immediate elongation of 
the (red) spring labeled G[. The (green) Voigt element 
then relaxes within a time scale of order r\\jG\. On very 
long time scales, the (blue) viscous element gives rise to a 
slow drift with velocity cr/773 (depicted by the finite slope 
of the blue line). 

In order to build a model that behaves like the Burger 
model under weak stresses but which presents a flow 
threshold like discussed earlier, one could combine the 
Burger model with the model presented in Fig.[T]for only 
one internal variable. We would thus obtain the model 
represented in Fig. [3b (which we already suggested as a 
generalization of our model [34[ ) . 

In the present work, we consider this combined model 
of Fig. |3Jd, but we focus on short and intermediate time 
scales where the highly viscous (blue) element has not 
moved yet. Because the blue element has not moved, 
it can be simply omitted. As for the two springs G\ 
and G[ and the viscous element 771, their behaviour can 
be reduced to that of a single spring in two limits. On 
time scales much shorter than r/i/Gi, the (green) Voigt 
element is blocked due to its viscous part. As a result, 
the whole system behaves like the red spring G[. At 
intermediate time scales (much longer than r\\jG\ but 
still with a blocked blue viscous element), the Voigt el- 
ement has relaxed. Hence, both springs are simply in 
series: they combine into a composite spring. In Fig. [3J;, 
we have represented such a model. The green and red 
spring represents either the red spring G\ (on short time 
scales) or the combination GiG[/(Gi + G[) (on interme- 
diate time scales). Meanwhile, the other creep elements 
provide both the stress threshold (solid friction 0-3) and 
the dependence on deformation rate (viscous element n' 3 ). 
We also include a general viscosity 772 as in the discussion 
of Paragraph III Al In this model, the spring and viscous 
elements must be understood as non-linear, unless stated 
otherwise. Apart from the viscous element 772, the model 
of Fig. [5J: is identical to the model that we constructed 
earlier and for which we had analysed the local, mean 
field behaviour 1341. 



A. General local rheological laws 

The relevant framework to describe elastic stresses in a 
flowing material is the Eulerian one, whether this mate- 
rial possesses clastical properties or not. Indeed, during 
the flow of a foam or an emulsion, even though clastic 
stresses exist, any reminiscence of a reference state dis- 
appears continuously due to plasticity. The Lagrangian 
description, which is based on maintaining the correspon- 
dance with such an intial state of reference, is formally 
equivalent, but less adapted conceptually and numeri- 
cally. 

Thus we attach the variables describing the material 
to a spatial grid (x, y, z), and they correspond to an in- 
stantaneous and local description in space. 

In this framwork, only two variables are relevant in a 
strictly mechanical context: the local velocity gradient 
Vt7(:e, y, z) and the local deformation state stored in the 
material j34[ (green-red spring in Fig.|3b), as described in 
continuum mechanics by the Finger tensor B(x, y, z) [49j . 
Note that when the material is at rest, B = I while 
the stored deformation, depicted schematically in Fig. [3J;, 
vanishes: e = 0. 

In this section, we describe our local rheological model. 
Note that in this local context, the global tensor V?j it- 
self has to be considered as an independent local three- 
dimensional tensorial variable, just as B, not as the spa- 
tial gradient of a velocity field. Only when we will turn 
to the description of a spatial system, see section IV El 
will the vector field v(x,y, z) be introduced. Meanwhile, 
tensors Vu and B will thus be the two variables of our 
local tensorial model. 

The elastic part of the stress, which goes throuth the 
spring in Fig.|3}: depends on the deformation according to 
the following relation, the most general one compatible 
with the symmetry constraints in three dimensions 3Jj: 



Ccl = ClQ I + Al B + Q2 B 



(3) 



where do, a\ and a-i are scalar functions of the invariants 
of the Finger tensor B. 

Turning to plasticity ((73 and n' 3 in Fig. [3]:), we only as- 
sume that every event of plastic relaxation is aligned with 
the stored deformation. The plastic creep should 
thus be similarly aligned. The most general form com- 
patible with the symmetry constraints is then: 



III. CONSTRUCTING THE TENSORIAL 
MODEL 



In the present Section, we will briefly recall how we 
built the rheological model [H[ of Fig. In particular, 
it is based on a general nonlinear description of elastic- 
ity and plasticity. Indeed, materials such as foams can 
locally undergo large elastic deformations — located far 
from the linear regime corresponding to small deforma- 
tions — before plastic flow occurs [U H|| . 



= b I + h B 



bo B 



(4) 



where bo, b\ and 62 are again scalar functions of the in- 
variants of the Finger tensor B. 

To complete the model, we gather together in a global 
viscosity term (which was noted 772 in Fig. all the dis- 
sipative phenomena which are present even in the absence 
of any plastic event in the foam. They occur for example 
at small scales: flows in films squeezed between bubbles 
or in Plateau borders. Wc simplify its description in se- 
lecting a Newtonian average viscosity rj s for these local 
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dissipative phenomena. The list of contributions to the 
stresses in the material is thus closed. We have: 



a = a,Q I + ai B + a 2 B 



(W+ v# 



(5) 



To take into account the incompressible character of 
foams and emulsions, we add an extra kinematic con- 
straint of strict volume conservation det(-B) = 1. Refer- 
ing to [HI for further details, we take it into account by 
using only the deviatoric part of the stress: 



dev(cr) 



a-Itr(cr). 



(6) 



The same constraint on plasticity gives the general 
form [H|: 

£>* = B-dev(f(B)) = h B-dev(B) + b 2 S-dcv(B 2 ), (7) 

where the scalar prefactors b\ and 62 are isotropic, and 
thus depend on the invariants of tensor B. 

In what follows, we will use a completely equivalent 
form of tensor which manifests more clearly that the 
dissipation is positive (see the discussion in |34|): 



D, 



A(B) 



B ■ Q(B) 



(8) 



where A{B) is a scalar isotropic function of B, r the 
characteristic time of the dissipative processes; moreover: 



5(B) 



dev [V(B) ■ dcv(cr c i)] 
tr[V(B) -devOd) • dev(cr 



(9) 



with V is a function of the form V{B) = b(B) B~ 2 + (1 - 
b(B)) B 2 [34|], where b is an isotropic function. In this 
expression, the total dissipation per unit volume is A(B) 
and can be chosen as positive. 

Eventually one gets the generic local rheological model: 



d_B 

lit 



Vu • B — B • Vv = -2D^, 



A(B) 



b-g(b), 



a = ao I + ai B + a-2 B 



Vv 



(10) 

(11) 
(12) 



where dB/dt = dB/dt+(v-V)B is the particulate deriva- 
tive of the Finger tensor. 



B. Complete spatial model 

As for any local rheological model, the previous equa- 
tions must be complemented by field equations which 
express force balance and mass conservation: 



dv 

P dt' Wp ' 



(13) 



dp 
Of 



+ V ■ (pv) = ^ +ptr^(Vv + Vv T ) = 0, (14) 



where / represents the external forces (per unit mass), 
and p is density. The incomprcssibility constraint gives 
here 



1 T 1 

V-iJ= tr-(Vu + Vv ) = 0. 



(15) 



As a result, the density p is simply transported by the 
flow: dp/dt = 0. In the remaining of this work, we fur- 
thermore assume that the density is homogeneous, hence 
it also remains constant: dp/dt = 0. 

Last assumption: we restrict ourselves to the Stokes 
regime, where inertial terms are all negligible in the mass 
conservation equation. Thus one obtains: 



-Vp. 



(16) 



The complete system of equations that we need to inte- 
grate numerically is thus: 



^-Vv-B-B-Vv T = -2D* 
dt p 

D* = 4&B.g{B), 

T 

1 T 

tr-(W+ Vv ) = 0, 
a = dev(tr) 
= dev (a 

V- a = -Vp, 



(17) 
(18) 
(19) 



dev { ai B + a 2 B 2 } + y(V« + Vw T ), (20) 

(21) 



The initial conditions that must be specified to solve 
the above system may merely consist in the values of 
tensor B over the entire sample. Indeed, the value of 
the velocity and pressure fields can be derived therefrom 
using equations (|21[) and (|20|) which, when combined, 
are similar to Stokes' equation, using the constraint of 
equation (jX9j) . 



C. Selection of a particular form of elasticity and 
plasticity 

1. Elasticity: Mooney-Rivlin model 

We have selected a usual form of incompressible elas- 
ticity which has been demonstrated to describe to a 
good approximation the nonlinear elastic behaviour of 
foams [50l IHlj : Mooney-Rivlin elasticity. The corre- 
sponding elastic energy per unit volume can be writ- 
ten [H: 

pE (B) = ^(l B -3) + ^-(n B -3) (22) 



where 



IIb 



tr(B) (23) 
J[tr 2 (B)-tr(B 2 )]=tr(B- 1 ) (24) 
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Going back the coeficients of Eq. (|3]), this corresponds to 
the following expressions: 



«i = ki + k 2 Ib 
a 2 = -k 2 



(25) 
(26) 



Following previous work refs. [50Ll5l|. we express the val- 
ues of k\ and k 2 using an elastic modulus G and an in- 
terpolation parameter a as follows: 



fci = aG 

k 2 = (l-a)G. 



(27) 
(28) 



In the foam modelling litterature, a value a = 1/7 is 
sometimes recommended (50l . [FH . Keeping in mind our 
perspective of discussing the conditions for the appear- 
ance of shear bands depending on parameter values, in 
sections IIVI and beyond, we prefer to keep the parame- 
ter a free, although we remain in the framework of the 
Moonlcy-Rivlin elasticity. 



2. Plasticity: yield stress fluid 

The particular form of plasticity explored in this work 
is based on a nonlinear threshold-like behaviour. Locally, 
the plastic reorganisation events only occur in the mate- 
rial when the stored elastic deformation reaches a critical 
value. We express this transition with a function W y (B) 
which vanishes linearly at the threshold: 



Wy{B) = 0, 



(29) 



with, in our case, W y (B) = pE(B) — K, where pE is the 
stored elastic energy per unit volume, and K a constant. 
In simple shear from a relaxed state, a y is the threshold 



stress: function W y vanishes. 



From the point of view of the plastic deformation rate 
tensor Dp 3 , we have the following expression [8j taking for 
A(B): 



A(B) = (pE(B) - K) Q{pE{B) - K), 



(30) 



where Q(x) = 1 when x > and Q(x) = elsewhere. 
We also set the following form for the polynom: 

V{B) = bB~ 2 + (1 - b)B 2 . (31) 

with b between and 1. Our final set of equations is thus: 

dB 



dt 



Vv-B-B-Vv T = -2 Dp 3 , 



(32) 



pE(B) - K 



Q(pE{B)-K)B-g(B), (33) 

7 

a = dev(cr) = dcv{(aG + (1 - a)Gtr(B)) B 
- (1 - a)GB 2 +7? S (W- 
V • a = -Vp, 

0. 



W T )/2} , 



1 t 
tr-(W+ Vv ) 



(34) 
(35) 
(36) 



3. Physical parameters and rheological model 

In order to be able to highlight physically relevant 
quantities, we use a non-dimensional form of the above 
system. The elastic modulus G is taken as the unit of 
stress, and the weak stress relaxation timescale r/ s /G as 
the unit of time, while B is already dimensionless: 



a = a/G, 
T = r, s t/G, 
B = B. 



(37) 
(38) 
(39) 



As a result, the various quantities are non- 
dimcnsionalizcd as follows: 



(40) 
(41) 
(42) 
(43) 
(44) 
(45) 
(46) 
(47) 



E(B) 


= PE(B)/G, 


K, 


= K/G, 


P 


= P/G, 


Mb) 


= A(B)/G, 


V(B) 


= V(B), 


G(B) 


- GG(B), 


w 


= (WG)W, 


D % 


= (r,./G)Dp 



The system of equations now reads: 



dB - - - - j 
— — = Vw ■ B + B ■ W 
dT 

tr(W+ W ) = 0, 



V • o 



a = CT e l 



<7 e i 



-Vp, 
w ■ 



Vv 



(a + (1 - a) tr(-B)) devB 
-(1 -a)dcvB 2 , 
VA(B)B-g{B), 
A(B) = (E(B) - K) Q(E(B) - 



D B 

p 



K), 



S(B) = |(I B -3) + i^(n fl -3), 



G(B)- 

V(B) 

* = 



dev [P(B) ■ gpi] 
tr [P(B) ■ a el ■ a e i)} 
(l-b)B 2 



bB- 2 



Gt 



(48) 

(49) 
(50) 

(51) 



(52) 
(53) 
(54) 

(55) 

(56) 

(57) 
(58) 



The new non-dimensional parameter "J defined in the 
last equation above reflects the ratio of the plastic flow 
rate (proportional to 1/t) to the viscoelastic flow rate 
(proportional to G/t] s ) when the other factors have the 
same order of magnitude. With the present choice for 
the magnitude of Dp 3 (with A proportional to the dis- 
tance from the threshold), that occurs when the stored 
deformation is typically twice the threshold deformation. 
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D. Simple shear flow 



In the remaining of this work, we address specifically 
the question of shear banding. For this purpose, we con- 
sider only simple shear flows. The velocity is oriented 
along axis x and varies along along axis y. The only 
non-zero component of the velocity gradient \7v is then 
dv x /dy. The entire system and flow are invariant along x 
and z. Besides, the force balance given by Eq. (|3"5j) then 



G, a 



implies that a xy and a yy are homogeneous at all times. 
In the axes x, y and z, the non-dimensionalized velocity 
gradient can thus be written as: 



Vv 



'0 T(y) 0^ 


,0 0; 



(59) 



where T(y) = { Vs /G)i(y). Let T = (r) s /G)i be the 
macroscopic value of the shear rate at the scale of the 
entire sample. We now have five non-dimensional pa- 
rameters: the plastic-to-viscoelastic flow rate ratio ty, the 
threshold /C, the Mooney-Rivlin parameter a, the param- 
eter b (which defines the tensorial form of the plasticity 
Q(B)), and the macroscopic shear rate T. 

For the sake of consistency, let us note that our earlier 
work [52j discussed parameters 



Q 



2 

r 



r 

Wc= - 
<3> 



(60) 
(61) 



instead of and T. 



E. Relation between model parameters and 
experimentally measurable quantities 

In fig. SJ we summarize the physical parameters in- 
cluded in our model. 

In experiments, easily accessible dimensional parame- 
ters are the viscosity r] s and the shear modulus G through 
linear rheology, as well as the threshold stress a y . More 
elaborate setups can yield the value of a. There are 
indications that a value a = 1/7 is relevant for liquid 
foams [HHll. 

Among our non-dimensional parameters, two can thus 
be determined easily: a and 1C. The latter is related 
to the energy at the threshold JC w ^(a y /G) 2 . As just 
mentionned T(y) = (j] s /G) 7(2/) is the normalised shear 
rate. Concerning ^ = i] s /(Gt) and b, no experiment to 
our knowledge is able to validate or invalidate the value 
of the plastic reorganisation time t at deformations be- 
yond the threshold, or the tensorial form of the plastic 
flow (here expressed in terms of parameter b). For the 
time being, we thus consider parameters and b as free 
parameters in any comparison of our model with actual 
data. 




FIG. 4: Simplified (scalar) picture of the main rheological pa- 
rameters. G represents the elastic modulus and a the relative 
weight of the tensorial components of the elastic deformation, 
see Eqs. (|27|) and (|28|) . The quantities a y and 1/t constitute 
a scalar representation of the creep defined by Dp 3 , and pa- 
rameter 6 is the equivalent of a for creep, see equation (|31l) . 
Finally, rj s is a viscosity that is independent of creep. 



IV. HOMOGENEOUS FLOW BEHAVIOUR IN 
LARGE AMPLITUDE OSCILLATORY 
EXPERIMENTS 

A. Method 

In the present section, we test the predictions of our 
model by comparing them to the most stringent available 
measurements in homogeneous flows, namely the large 
amplitude oscillatory experiments conducted recently by 
Rouyer et al. [11]. We have simulated oscillatory shear 
flow with the local model (no dependence on coordinate 
y, i.e., homogeneous flow). In other words, in Eq. (|59p . 
we choose 

t(y,t) = t(t) = -ojTq cos(wi), (62) 
which corresponds to the oscillating shear deformation: 
r(y,t) = r(t) = r sm{ujt), (63) 

B. Typical behaviours 

Figs. O and [5] show, in the form of Lissajous curves, 
the full response of the present model in the (amplitude, 
frequency)-plane. The normalized stress and strain re- 
sponses make it easy to discriminate between plastic, 
elastic, or viscous behaviours of the model. Let us ra- 
tionalize them in terms of the scalar diagram of Fig. |3J: 
discussed above in Section III Bl 

A pure elastic behaviour corresponds to an ellipse 
squeezed into a straight line spanning the diagonal of the 
diagram. This is obtained at low frequencies and ampli- 
tudes. Indeed, deformation rates are then small at all 
times, hence the viscous elements in diagram [3J: play no 
role. Meanwhile, because deformations remain small, the 
threshold of the solid friction clement is never reached. 
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■Q 


o 

(y 

/ 

/ 


o 

u 


o 

o 
o 


o 
o- 

o 


0.01 0.1 1 

UI 

Shape of normalized Lissajous c 


10 100 

urves (imposed strain) 


o 

Z7 
/ 

0.01 


o 
/ 

0.018 


o 

/ 

/ 

0.032 


o 

/ 

/ 

0.056 


/ 

/ 

0.1 



FIG. 5: Large amplitude oscillatory simulations: shape of F(t) 
versus cr xy (t) (Lissajous) curves, obtained for a = b = 1/7, 
* = 0.1 and K. = 1. Top: range from ui = 0.01 to cj = 100 
and To = 0.1 to To = 100. Bottom: zoom on a more restricted 
range of parameters. 



As a result, the spring alone provides the mechanical re- 
sponse of the system. 

For the same reason, an elasto-plastic behaviour is ex- 
pected at low frequency yet large amplitude, since the 
stress threshold is then reached. A purely elasto-plastic 
behaviour, as predicted by a scalar model, would corre- 
spond to a sharp-cornered parallelogram with two hori- 
zontal sides corresponding to the yield stress. The results 
of our simulation at low frequency and large amplitude 
differ from this simple picture in the same way as experi- 
mental data by Rouyer et al. [38j , namely with two main 
features: (i) the "plastic part" of the Lissajous curve ex- 
hibits a slightly negative slope, and (ii) the transition to 
plasticity is progressive rather than sharp (blunt corners) . 
Feature (i) corresponds to the weakening of the viscous 
component when the deformation rate decreases along 
the sinusoidal applied deformation. As for the latter fea- 
ture, it can result either from viscosity being combined 
with plasticity (as in the present model [34]) or from a 



ui = o.oi, r = l to ioo 




-1 -0.5 0.5 1 



r(t)/r„ 



ui = 0.01, To = 1.4 and 1.8 




-2 -1.5 -1 -0.5 0.5 1 1.5 2 



r(t) 

FIG. 6: Large amplitude oscillatory simulations: shape of 
0xy(t) versus T(t) (Lissajous) curves, obtained for a — b = 
1/7, ui = 0.01, $ = 0.1 and K, = 1. Top: transition from elas- 
tic to plastic behaviour as the amplitude is increased. At very 
large amplitudes, an overshoot is apparent like in continuous 
shear situations. Bottom: in a slightly plastic situation (am- 
plitude 1.8), it takes several cycles before the system behaves 
in a periodic manner. 

progressive onset of plasticity [3(| . 

As compared to the results by Rouyer et al. [3a |. our 
model additionally exhibits an overshoot at very low fre- 
quency and large amplitude. Because such a regime is 
very similar to slow, continuous shear, this response can 
be understood [35| as a tensorial effect combining the 
saturation arising from plasticity and the rotation con- 
tained in shear (this point is further discussed at the end 
of Section EH). 

This transition between a purely elastic response at low 
amplitude and an elasto-plastic response at higher am- 
plitude is best illustrated by the top part of Fig. [5] The 
curves are normalized for clarity, but the actual maxi- 
mum slope in each curve is essentially identical and is 
given by the shear modulus G. By contrast, the value of 
the stress in the most horizontal regions of the curve re- 
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time / period 

FIG. 7: Comparison of model (curves) with experiments 
(points). Normalized a(t) curves for six values of the am- 
plitude To ranging from 0.055 to 1.2, shifted vertically for 
clarity. Parameter values are a = 0.14, K, = 0.04, b = 0.14, 
$ = 27. For each value of the amplitude F we had to select 
a different value for the frequency u> and for the modulus G. 
The (F, uj, G) values are: (0.055, 0.2, 282), (0.15, 0.4, 246), 
(0.25, 0.9, 186), (0.43, 0.5, 154), (0.72, 0.5, 146), (1.2, 0.4, 
149). 

fleet both the solid friction element and both viscous ele- 
ments in diagram^. The bottom part of Fig. [5] presents 
results slightly below and slightly above the plasticity 
threshold. It shows that weak plasticity causes the de- 
formation cycle to slowly drift towards a limit cycle that 
differs from the elastic cycle. This very drift, when con- 
tinued along a longer cycle, is in fact at the origin of the 
overshoot mentioned above. 

If we now turn to higher frequencies, the deformation 
rate becomes large. As a result, viscous elements now 
play a role and may even become dominant. Correspond- 
ingly, the Lissajous curve tend to become an ellipse whose 
axes lie along those of the figure. That is particularly 
clear on diagram [3J: at high frequency with a large am- 
plitude, but the trend is very obvious at high frequency 
and low amplitude, and is also discernable at low fre- 
quency and large amplitude. 



C. Comparison with experiments 

The data obtained by Rouyer et al. [38| correspond to 
a fixed frequency and different applied strain amplitudes 
(from 0.055 to 1.2). We have integrated [59[ the equa- 
tions of the present model with the same amplitudes and 
plotted them together with data. The results are pre- 
sented in figure [7] (for clarity, stress curves are presented 
normalized) . 

Note that in order to obtain a reasonable agreement 
of our model with the data, we had to artifically choose 
different values of ui and G for each strain amplitude, 



while k and iff could be kept constant. This unsatisfying 
ad hoc parameter adjustment shows the limits of this 
model in describing the behaviour of the foams studied 
by Rouyer et al. 

V. SHEAR BANDING STUDY 




FIG. 8: Typical form of a stationnary flow curve giving the 
dependence of the shear stress on the local shear rate. a y 
is the yield stress as measured under imposed stress, and j c 
the corresponding shear rate. A macroscopic shear rate j A 
smaller than will not necessarily lead to a homogeneous 
velocity profile P A , with the expected stress a A : the flow can 
separate into a blocked region and a flowing region (profiles 
P A or P A ). The local shear rate is then faster (7^2 > j A and 
7A3 > 7 A )i which corresponds to a higher stress (a A > a A 
and a A > a A ). Besides, for an average shear rate 7 s greater 
than 7c, the flow is homogeneous again, which corresponds to 
the expected stress a B (greater than a y ). 



A. Stationary flow curve and inhomogeneous flow 

Let us now turn back to shear banding. For such a pe- 
culiar flow to be observed, the same material submitted 
to the same shear stress o~ xy must be simultaneously in 
two different deformation states. As discussed for many 
years for various complex fluids 0, IBS , , a mathemat- 
ical condition for this to be possible is the existence of 
an unstable zone in the local flow curve (700(7) of the 
material: it must be non-monotous. 

In the case of foams, nevertheless, such an unstable 
portion in the flow curve itself does not exist: how can 
shear bands with different shear rates coexist? 

Foams and emulsions arc instances of yield stress flu- 
ids, so that there exists a minimal value o~ y of the stress 
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a xy below which no stationary flow occurs. Now when 
we shear the material, imposing the shear rate, the ma- 
terial has to flow, even for very small 7. The intrinsic 
flow curve thus possesses an extrapolation in stress when 
7 — > 0. Let's denote it by ad- Note that a y and ad 
pertain to the local rheology curve 1700(7), n °t to the ef- 
fective, macroscopic stationary curve as can be measured 
for example in a rheomctcr. In this discussion the flow is 
homogeneous. But the relative values of ad and a y , per- 
taining to the local flow curve, will give us hints about 
possible conditions for shear banding. 

Let us now consider the result of a measurement made 
on a sample of this material, sheared in a (parallel or very 
low curvature) Couette cell under imposed shear rate. If 
> a y , all parts of the sample will flow, even at low 
shear rates, since the corresponding stress is necessarily 
everywhere greater than the yield stress. As mentionned 
before, since the flow curve has no intrinsic instability 
for higher 7 values, no mechanism is available for shear 
banding. 

The situation is different if ad < a y . If we put on the 
same graph the yield stress a y and the intrinsic station- 
ary flow curve (figure [5]), it is immediately apparent that 
this configuration allows for the coexistence of zones un- 
dergoing shear at rates 7 such that < 7 < j c , and of 
blocked zones remaining in the elastic regime at 7 = 0. 
The mechanism is essentially the same as in the classi- 
cal case of instability in the flow curve (see figure [8|). Of 
course, as soon as 7 > j c all regions flow, since 7 > j c im- 
plies that some regions flow faster than j c . The stress in 
these regions, as given by the flow curve, has to be above 
the yield stress a y . And since the shear stress is the same 
in the entire material, all regions support a stress greater 
than a y and no region can be blocked. 



Large elastic deformations versus extra dynamic 
variables 



But is the situation where ad 7^ a y actually possible? 
In various complex fluids, the answer is known to be yes. 
The usual explanation of such a flow curve is to invoque 
an internal extra variable (of a structural nature in gen- 
eral) which is coupled to the flow. As an example, in 
a simplified vision, this extra parameter can take one of 
two values: flowing or non-flowing. Thus the stationnary 
curve extrapolating to ad at low 7, and the yield stress 
value a y correspond in reality to two different materials, 
hence ad and a y can differ. 

But as mentioned in Section II B( foams differ from 
many other complex fluids in that the deformation that 
must be reached to trigger plastic flow is large. This fea- 
ture turns them into an intrinsically tensorial material. 

In a stationary situation where shear banding is 
present, stress conservation implies that the shear stress 
a xy is constant along the direction of the velocity gra- 
dient, as well as the stress component a yy . By contrast, 
the extra components of the stress a xx (y) and a zz (y) may 



vary in an arbitrary manner along the direction of the ve- 
locity gradient. Among these, a xx (y) is present even in a 
purely 2D system. These extra components will qualita- 
tively play the same role as an extra structural variable 
in changing the local nature of the material, when viewed 
as a ID material (along direction y). 

But that only explains how it is possible for shear bands 
to be present. The reason why the flow curve actually ex- 
trapolates below the yield stress at vanishing shear rates 
{ad < a y ) in some tensorial models, thus allowing shear 



banding, has been shown by Raufastc et al. [35[ : as long 
as the material remains elastic, the local deformation ten- 
sor is transported by the shear flow along a path that is 
not locally aligned with itself: the principal axes of the 
particulate time-derivative of the deformation do not co- 
incide with those of the deformation. Hence, once plastic- 
ity is triggered, it alters the deformation evolution until 
it progressively reaches the locus where it is aligned with 
its transport under shear. At least in simple examples of 
elasticity and plasticity, this migration from the clastic 
path to the asymptotic locus is the origin of the shear 
stress overshoot observed during transients [35|. When 
the plastic flow is triggered rather abruptly, the system 
is still elastic just before the maximum of this overshoot, 
and the asymptotic shear stress value can then lie be- 
low the last elastic shear stress value. In other words, 
ad < a y . 



C. History dependent shear bands 

Despite some similarities, the analogy with systems 
characterised by unstable flow curves has some limita- 
tions. In the case of yield stress fluids, there is no un- 
stable range in 7, which would impose phase separation 
between two phases at different flow rates. Shear bands 
are possible but not necessary. Also, no lever rule-like 
criterion can exist to select the relative fraction of the 
different bands, as have been argued in some fluid sys- 
tems 0,[5| . 

Rather, the initial distribution of a xx {y) and a zz {y) 
in the material will be of primary importance in the ap- 
pearance of shear bands even though the material per se 
remains homogeneous. In other words, it is the material 
history that will lead to a particular flow profile. We will 
see that the initial distribution of stress in the material 
will determine the band structure. 



D. OD flow curve and shear banding criteria 

As long as the flow in the material is homogeneous, 
a local rheological model will be sufficient to describe it. 
We begin by showing the typical flow curve corresponding 
to our model (fig. [9]). Note that this flow curve is obtained 
under applied shear rate conditions. 

As can be observed, the conditions described in the 
introduction for the appearance of shear bands are 
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applied to the system in the stationary state would pre- 
cisely correspond to the plastic threshold: cr^2 at (7c) = a y 
Such a procedure is natural, but requires successive sim- 
ulations of the system for a large number of 7 values. 

We have used a more direct approach [3J] to obtain 
ad and 7 C (see Appendix). This method relies on the 
description of the evolution of the system in terms of 
independent eigenvalues pi and of tensor B (see fig- 
ure [TDJ). 

With the help of this procedure, the three observables 
which are important for the prediction of shear bands, 
a y , <Jd, and 7c, are obtained directly without the need for 
simulating separately all the points along the stationary 
flow curve. 
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FIG. 9: Top: typical stationary flow curve. Points corre- 
sponding to Gd, <J y and 7 C are reported on the curve. Bottom: 
stress time evolution for different imposed shear rates below 
or above 



FIG. 10: Form of the stored elastic deformation in the course 
of an experiment and in the stationary regime, for three dif- 
ferent values of the shear rate 7. The axes are the first two 
eigenvalues, (3\ and P2, of tensor B. 



fulfilled: stress a c i is smaller than the static yield stress 
tjy . In the shear rate range between and j c , the system 
has the possibility to split the average shear rate 7 in 
different proportions of blocked and flowing bands. 

Thus, in the homogeneous case, for any values of the 
parameters JC, a, b and T, we can use the local rheo- 
logical model to calculate the static and dynamic thresh- 
olds, a y and o"d, and the critical shear rate 7 C . Following 
the line of reasoning developped in the introduction, we 
can then predict the range of imposed shear rates [0, j c ] 
inside which shear bands are possible. 

The value of a y can be obtained easily by simulating 
the system in the elastic regime (Dp 3 = 0) up to the 
threshold (W y (B) — 0), which corresponds to a state 
of the system characterized by eigenvalues /3f and f3\ of 
tensor B, a state for which a v can be calculated. 

The different stationary state values of the shear stress 
could then obtained independently by continuing the sim- 
ulation beyond the threshold in the plastic regime for 
each value of 7, waiting for the stationnary value of the 
system (dB/dt ps 0). The dynamic threshold a 4 would 
then correspond to the limit of <J\2 for small 7. The 
critical shear rate 7 C would be obtained when the stress 



E. Spatial (ID) simulations of stationary flow 
regimes 

As already mentioned, the discussion in the previous 
paragraph only provides necessary conditions for the ap- 
pearance of shear bands. In the ID simulations that will 
be discussed in the present section, flow inhomogeneities 
will indeed sometimes emerge in a full spatial simulation 
of our tensorial model in ID spatial dimension plus time. 

We simulate the full tensorial model in ID using the 
equations of paragraph IHI C 31 The technical details of 
the numerical scheme can be found in [52{ . From a 
numerical point of view, let us remark in particular that 
we have checked the the grid used in the discretisation of 
the equations is fine enough for all simulations presented 
here. 



1. Discussion of the conditions for inhomogeneous flow 

The model that we simulate only contains material pa- 
rameters that are homogeneous in the sample. Hence, if 
the initial conditions of the flow are also homogeneous, 
the entire evolution will remain homogeneous. Although 
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performing a ID simulation as a set of partial differential 
equations, we would obtain the exact same results as in 
the previous Section. 

In other words, since the parameters of the model do 
not vary in space, shear bands can only appear if initial 
conditions are, in one way or another, inhomogeneous. 

Of course, as mentioned in the Introduction, inhomo- 
geneities could appear in a natural way through an extra 
state variable coupled to the flow, such as the concen- 
tration. This variable could then vary in space and be 
coupled to the flow. Concerning concentration (a con- 
served variable), let us mention dilatancy phenomena, 
imagined for foam [13, HH , observed experimentally 
and interpreted in a geometrical manner [2^, [57| ■ Aligne- 
ment (a non-conserved variable) is another possibility. It 
has been invoked in the case of wormlike micelle or rigid 
rod solutions (55| . 

Here, we focus on inhomogeneous static strain/stress 
initial conditions, without invoking additional variables, 
and we will show that they can induce the appearance of 
persistent inhomogeneities in the flow profile. 

The reason for which these initial strain inhomo- 
geneities can induce the appearance of blocked bands 
can be qualitatively undestood by considering the flow 
threshold JC. Indeed, the stresses generated by the shear 
combine with the initial stress distribution due to strain 
inhomogeneities. Depending on its orientation, the ini- 
tial stress thus precipitates or delays the triggering of the 
plastic flow. 



2. Initial inhomogeneous strain distribution 

First, the existence of stress inhomogeneities stored in 
the system before it is set into motion is physically well 
motivated. For instance, introducing a foam sample into 
an apparatus requires non-homogeneous flows. Inhomo- 
geneous stresses will build up in the sample very likely, 
except if particular care is taken, such as a slow, in situ 
drying of an initially wet foam. 

We will always assume that the initial state is at rest, 
that is, that the elastic stresses are at equilibrium in the 
sample. However, even when this equilibrium is imposed, 
there exist a large set of possible initial spatial distribu- 
tions of stresses and strains. For example, if the system 
is invariant in the xz plane of the shearing walls, some 
components of the stress must be homogeneous. That is 
the case for a xy , a yy and a yz . The other stress compo- 
nents, however, can freely vary as a function of y as long 
as they remain constant in each xz plane. It thus cor- 
responds to a ID inhomogeneity in the direction of the 
velocity gradient. 

In this paragraph, we examine a very simple case of 
initial condition, with uniaxial extension along axis x, 
with B xx = B xx (y) and B yy = B zz = l/y/B xx . We used 
a simple monotonic function: 

B xx = l.l + eyP(l-(l-yf). (64) 



In practice, in order to prepare a sample in such a 
state, one must compress the foam in a non-homogeneous 
manner. Typically, a block of foam with a trapezoidal 
shape forced to take a rectangular shape will undergo 
this kind of strain inhomogeneity. In this context, we 
cannot comment on any relation between these strain 
inhomogeneities in our continuum model and the local 
structural disorder existing at the bubble level, as this 
disorder is averaged out in our continuum description. 
In particular, there is no clear structural interpretation 
of the amplitude e of the strain inhomogeneities in the 
prepared sample. 



3. Characterizing the inhomogeneous flows 

A typical sequence of velocity profiles obtained in our 
numerical simulations displays as follows. The velocity 
profile is initially homogeneous. It remains homogeneous 
as long as the entire sample is in the elastic regime. The 
regions where the initial stress is the highest in the direc- 
tion of the applied deformation reach the threshold first. 
The average shear rate being constant, this onset of creep 
leads both to a higher shear rate in the creeping regions 
and to a lower one in the others. The high shear rate then 
induces the saturation of the stress due to creep, and the 
shear becomes blocked in the region below the threshold. 
In the stationary regime, a blocked band coexists with a 
sheared band at the same shear stress. 

In the corresponding transient regime, non-trivial phe- 
nomena may appear, especially at the boundary of the 
blocked zone. Transient negative local shear rates are 
observed due to stored elastic stresses. 

Let us now address the characteristics of the stationary 
velocity profile, again from the behaviour of the local 
rheological model. 

The first feature of interest is that in the flowing re- 
gions, the velocity profile is linear, that is, the shear rate 
is uniform. That can be understood in the following 
manner. All the regions which, in the stationary state, 
respond through a non-zero shear rate, correspond to a 
point located on the stationary flow curve in the f5\-$i di- 
agram of figure 1101 Each point of this curve corresponds 
to a different shear stress. Thus, since each layer of the 
flow undergoes the same shear stress, they all actually 
correspond to the same point on the curve and thus re- 
spond through the same shear rate. 

A second feature results from the fact that in space, 
while o~ xy and a yy are continuous, o~ xx and a zz can per- 
fectly be discontinuous. That is precisely the case at the 
boundary between a shear and a blocked region. This is 
the flow counterpart of the discontinuity in the /3i-/32 di- 
agram, between the points below the threshold and the 
point with a stationary shear rate that corresponds to 
the flowing region. Actually, the only coupling between 
the different layers comes from the fact that (i) o~ xy and 
o~ yy must be every where the same, and (ii) the integral 
of 7 over the gap thickness is fixed by the imposed wall 
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velocity. As a consequence, in an inhomogcncous flow, 
the organisation of the blocked and flowing layers is not 
unique: any permutation of the layers is actually possi- 
ble. Again, initial conditions decide upon the particular 
structure adopted by the flow. Two initial conditions cor- 
responding to permutated layers would lead to the same 
permutation in the stationary flow structure. 



Parameters affecting the existence of blocked 
bands 



In this section, we want to describe, within the pa- 
rameter space (iff, JC, a, b, T) the regions inside which 
shear bands are possible. These domains will be repre- 
sented through sections in five different planes: (JC,T), 
('J/,]?), (a, b) and (iff, JC). The results are presented in 
figures HU [H [HI [H and HU 

As will be discussed below (paragraph IV F 6[) . the 
choice of the initial conditions can have a crucial impact 
on the existence of shear bands. A complete investiga- 
tion of the model would therefore require a very thor- 
ough exploration not only of the parameters (iff, JC, a, 
b, r) but also of the shape and amplitude of the initial 
strain profile. In order to favour the appearance of shear 
bands without needing to refine the exploration of var- 
ious strain profiles, we selected a very large amplitude 
e = 500% for the shape mentioned in Eq. (fM)) . This 
applies to Figs [TT1fT5l 



1. (JC, if 1 ) plane 

In the (JC, iff) plane, bands are predicted by the local 
model for large values of iff and JC. Note that in the figure, 
zones where no bands can appear are denoted by blue 
triangles. Small values of iff correspond to a situation 
where the relaxation time in the absence of plasticity 
is far smaller than the relaxation time r corresponding to 
the plasticity. It is thus a regime dominated by the fluid 
viscosity, where the creep plays no role. As iff increases, 
the creep becomes dominant, and shear bands can appear 
for lower values of the threshold (JC) (figure fTTj) . 

As expected, the shear bands observed in the simu- 
lations appear only in regions authorized by the scalar 
model. 

The green dots correspond to values for wich the scalar 
model allows the presence of shear bands, wich are not 
observed in the simulations for a specific set of initial 
conditions. As will be commented further on, the extent 
of this green zone depends on these initial conditions, 
demonstrating one of the main points of this work. 



2. (JC, V) plane 

In the (JC, T) plane, bands should be predicted for 
small values of T (due to the small velocities wich explore 



0.05, a = 6= 1/7 
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FIG. 11: Comparison of 0D and ID simulations in plane (JC, 
iff), using a = b = 1/7 and T = 0.05. Blue triangles indicate 
values for which the 0D model allows only uniform flow. Red 
squares indicate values for which shear banding was obtained 
in the ID simulation for the initial conditions chosen. Green 
disks indicate additional values for which the 0D model allows 
shear banding. 



regions of the flow curve close to the origin in T), and for 
large values of JC. Indeed, in that case the static threshold 
is large which favours bands since they are possible below 
this threshold (figure IP2j) . 

Concerning the relation between the values predicted 
for shear bands in the OD model and the observations 
in the ID simulations, the same remarks hold as for the 
previous paragraph. 



* = 0.1, a = b= 1/7 
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10 
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FIG. 12: Comparison of 0D and ID simulations in plane (r, 
JC), using a = b = 1/7 and ^ff = 0.1. Blue triangles indicate 
values for which the 0D model allows only uniform flow. Red 
squares indicate values for which shear banding was obtained 
in the ID simulation for the initial conditions chosen. Green 
disks indicate additional values for which the 0D model allows 
shear banding. 
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3. (*,r) plane 



* = 0.1, r = 0.05, b=l-a 



Observations in this plane corroborate the analysis in 
the two previous planes : bands are allowed (and are 
observed) for low T values and high VP values (figure IT31 . 
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FIG. 13: Comparison of 0D and ID simulations in plane (f , 
9), using a = b = 1/7 and JC = 1.0. Blue triangles indicate 
values for which the 0D model allows only uniform flow. Red 
squares indicate values for which shear banding was obtained 
in the ID simulation for the initial conditions chosen. Green 
disks indicate additional values for which the 0D model allows 
shear banding. 



4- (JC, a) plane 

Again bands appear for large values of JC. The in- 
fluence of the a parameter is far more subtle to assess, 
being related to non-trivial tensorial effects of the elastic 
(a) and plastic (b taken as 1 — a here) terms. 

The same remarks hold concerning the correlation be- 
tween the 0D model and ID simulations. 



5. (a, b) plane 

Finally, in the (a, b) plane, one is again confronted with 
3D effects which are difficult to discuss in intuitive terms 
(figure I15[) . The way the elasticity (parameter a) and 
the plastic deformation rate (parameter b) are coupled 
in a tensorial way affects the critical rate 7 C and can be 
enough to eliminate all possibilities of shear bands. 



6. Dependance on the initial conditions 
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FIG. 14: Comparison of 0D and ID simulations in plane (JC, 
a), using b = 1 — a, T = 0.05 and $ = 0.1. Blue triangles 
indicate values for which the 0D model allows only uniform 
flow. Red squares indicate values for which shear banding 
was obtained in the ID simulation for the initial conditions 
chosen. Green disks indicate additional values for which the 
0D model allows shear banding. 



However, the boundary of these regions do not coin- 
cide. In fact, for the same values of the parameters, the 
extent of the banding zone depends crucially on the initial 
conditions, while always remaining in the region allowed 
by the rheological model. In other words, the behaviour 
of the system is history dependent, a feature realized in- 
dependently in a recent work on a related tensorial model 
with plasticity fo8l |. 

To illustrate that, we have varied both the form and 
the amplitude of the spatial modulation of the initial de- 
formation. In Figures [TT1I151 the initial profile was given 
by the non-linear form of Eq. ([64)) with e = 500%. By 
contrast, in Figure 1161 for the upper graph we chose a 
simple sigmoidal profile centered around a selected alti- 
tude y , 



B, 



1 



ey 



Vo 



1 



y o 



(65) 



and for the lower graph we chose a step-like function, 
both with e = 10%. Comparing figures [T21 and [TBI shows 
that the parameter domain where shear bands actually 
appear can depend in a non-trivial manner not only on 
the shape but also on the amplitude of the initial strain 
profile. 

Actually, we expect that a thorough exploration of the 
region authorized for the shear bands could be achieved 
through a very fine adjustment of the initial condition 
profile for each set of parameters. 



We have always observed that the regions in which 
blocked bands actually appeared in the ID simulations 
are strictly included, as expected, in the regions autho- 
rized by the local rheological model. 



VI. CONCLUSION 

The present study on shear bands in liquid foams was 
conducted on a rheological model (3~H whose predictions 
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FIG. 15: Comparison of 0D (top) and ID (bottom) simula- 
tions in plane (a, b), using T = 0.05 and = 0.1. Following 
the indications of Fig. [14] we chose K, = 0.3 for 0D simula- 
tions and /C = 0.8 for ID simulations. The second diagonals 
(b — 1 — a) in the present diagrams correspond to vertical lines 
in Fig. 1141 Blue triangles indicate values for which the 0D 
model allows only uniform flow. Red squares indicate values 
for which shear banding was obtained in the ID simulation for 
the initial conditions chosen. Green disks indicate additional 
values for which the 0D model allows shear banding. 
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FIG. 16: Initial conditions dependency in the plane (f , JC), 
using a = b — 1/7 and 'J? = 0.1. In the upper graph we 
considered sigmoidal initial conditions, and in the lower graph 
step-like ones. In both cases, the amplitude of the strain 
inhomogeneities was reduced to 10% in amplitude, whereas 
Fig. 1121 corresponded to 500% to enhance the effect. 



we here compare to existing rheological measurements 
under large amplitude oscillations (see Section HV CI) . 

The shear bands obtained with the model in parallel 
Couette geometry display several somewhat unusual fea- 
tures. 

1. The shear bands depend on the initial conditions. 
More generally, the stationary state is history de- 
pendent: only the flowing regions coincide with the 
simple linear velocity profile obtained in the case of 
a stationary homogeneous flow. 

2. The response of the model in a stationary homo- 
geneous flow is continuous when approaching zero 
shear rate, with no forbidden region below some 
finite shear rate. 

3. The model does not contain any non-conserved 



(structural) order parameter. 



This study thus shows that shear bands can arise nat- 
urally in a fully tensorial rheological model. This departs 
from most works in the shear banding community which 
put less emphasis on the tensorial character of the vari- 
ous models. Here, the appearance and persistence of the 
bands result from the combination of the initial condi- 
tions and the difference between the static and the dy- 
namic flow thresholds (in shear geometry), which itself 
arises from the tensorial character of the model. The 
shear rate is discontinuous at the boundary between the 
flowing and blocked regions, but the value of the shear 
rate near the boundary as well as the band widths depend 
on the sample history and preparation. 
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Appendix A: Direct method for obtaining the 
stationary state in the local rheological model 

Let us start from the point (/3f , and follow the 
plasticity threshold W y (B) = until the stationarity con- 
dition is fulfilled. This condition can be expressed using 
the following observation: in the stationary reg ime, there 
is no plastic flow in the vorticity direction [34| . In other 
words, the third eigenvalue of tensor G(B) is zero: 



g 3 {fa,fa) = Qz{fa,fa,fa) = Q, 



(Al) 



with fa = We thus directly obtain the dynamic 

threshold (/3f, (3^) of the system. We then follow the 
same stationarity condition g 3 (fa, fa) = until we reach 
the desired shear stress cr 12 = cr y . We thus directly obtain 



the stationary state (/Jf 3 , fi^) that corresponds to the 
critical shear rate j c . In practice, we follow the threshold 
curve using W v (fa,fa) = W y (fa, fa, fa) = (with fa = 
p\ ) by integrating the following differential system: 



dfa 


dWy 


dWy 


1 


dWy 


dt 


dfa 


dfa 




dfa 


dfa 


dWy 




1 


dWy 


dt 


dfa 


dfa 




dfa 



(A2) 
(A3) 



where the sign of Ew 



follow the curve W y in the desired direction 



± 1 is chosen in such a way as to 
Similarly, 

we follow the curve of stationary states, g 3 {fa,fa) = 
G 3 (fa,fa, fa) = by integrating the following differential 
system: 

aq. Pi„ Fir. 1 fir 

(A4) 
(A5) 



dfa 


dg 3 


dG 3 
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dg 3 


dt 


' dfa 


dfa 


/3i/? 2 2 


dfa 


dfa 


dg 3 


dG 3 


1 


dQ 3 


dt 


' dfa 


dfa 




dfa 



where the sign of e ga = ±1 is chosen in such a way as to 
follow the curve 173 = in the desired direction. 
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